Load Data
library(gdata)
dataset <- read.delim("raw_data/Figure3D.txt", stringsAsFactors = FALSE)
dataset$genotype <- factor(gsub(" ","",dataset$genotype))
dataset$experiment <- factor(rep(paste0("exp", 1:(length(dataset$genotype)/length(levels(dataset$genotype)))), each= length(levels(dataset$genotype))))
dataset$UID <- factor(paste(dataset$experiment, dataset$genotype))
kable(dataset, row.names = F)
| WT |
8.366667 |
35.60000 |
29.56667 |
21.23333 |
exp1 |
exp1 WT |
| PARP1KO |
6.700000 |
57.23333 |
46.26667 |
45.96667 |
exp1 |
exp1 PARP1KO |
| ALC1KO |
6.533333 |
56.96667 |
59.76667 |
55.66667 |
exp1 |
exp1 ALC1KO |
| ALC1KOPARP1KO |
6.400000 |
55.26667 |
43.00000 |
46.40000 |
exp1 |
exp1 ALC1KOPARP1KO |
| WT |
7.300000 |
36.13333 |
29.73333 |
23.00000 |
exp2 |
exp2 WT |
| PARP1KO |
7.966667 |
48.96667 |
49.60000 |
43.93333 |
exp2 |
exp2 PARP1KO |
| ALC1KO |
10.600000 |
54.16667 |
59.83333 |
57.93333 |
exp2 |
exp2 ALC1KO |
| ALC1KOPARP1KO |
8.300000 |
53.43333 |
46.23333 |
42.43333 |
exp2 |
exp2 ALC1KOPARP1KO |
| WT |
7.233333 |
34.06667 |
35.43333 |
20.66667 |
exp3 |
exp3 WT |
| PARP1KO |
7.233333 |
39.76667 |
41.26667 |
33.56667 |
exp3 |
exp3 PARP1KO |
| ALC1KO |
7.233333 |
54.46667 |
52.10000 |
55.13333 |
exp3 |
exp3 ALC1KO |
| ALC1KOPARP1KO |
6.433333 |
42.43333 |
42.90000 |
35.96667 |
exp3 |
exp3 ALC1KOPARP1KO |
library(reshape2)
# reshape to long format
dataset <- melt(dataset, variable.name = "treatment", value.name = "Outcome")
dataset$genotype <- relevel(dataset$genotype, ref = "WT")
dataset$UID <- relevel(dataset$UID, ref = "exp1 WT")
dataset$time <- gsub("NT","0",dataset$treatment)
dataset$time <- gsub("X|h","",dataset$time)
dataset$time <- as.integer(dataset$time)
kable(dataset, row.names = F)
| WT |
exp1 |
exp1 WT |
NT |
8.366667 |
0 |
| PARP1KO |
exp1 |
exp1 PARP1KO |
NT |
6.700000 |
0 |
| ALC1KO |
exp1 |
exp1 ALC1KO |
NT |
6.533333 |
0 |
| ALC1KOPARP1KO |
exp1 |
exp1 ALC1KOPARP1KO |
NT |
6.400000 |
0 |
| WT |
exp2 |
exp2 WT |
NT |
7.300000 |
0 |
| PARP1KO |
exp2 |
exp2 PARP1KO |
NT |
7.966667 |
0 |
| ALC1KO |
exp2 |
exp2 ALC1KO |
NT |
10.600000 |
0 |
| ALC1KOPARP1KO |
exp2 |
exp2 ALC1KOPARP1KO |
NT |
8.300000 |
0 |
| WT |
exp3 |
exp3 WT |
NT |
7.233333 |
0 |
| PARP1KO |
exp3 |
exp3 PARP1KO |
NT |
7.233333 |
0 |
| ALC1KO |
exp3 |
exp3 ALC1KO |
NT |
7.233333 |
0 |
| ALC1KOPARP1KO |
exp3 |
exp3 ALC1KOPARP1KO |
NT |
6.433333 |
0 |
| WT |
exp1 |
exp1 WT |
X2h |
35.600000 |
2 |
| PARP1KO |
exp1 |
exp1 PARP1KO |
X2h |
57.233333 |
2 |
| ALC1KO |
exp1 |
exp1 ALC1KO |
X2h |
56.966667 |
2 |
| ALC1KOPARP1KO |
exp1 |
exp1 ALC1KOPARP1KO |
X2h |
55.266667 |
2 |
| WT |
exp2 |
exp2 WT |
X2h |
36.133333 |
2 |
| PARP1KO |
exp2 |
exp2 PARP1KO |
X2h |
48.966667 |
2 |
| ALC1KO |
exp2 |
exp2 ALC1KO |
X2h |
54.166667 |
2 |
| ALC1KOPARP1KO |
exp2 |
exp2 ALC1KOPARP1KO |
X2h |
53.433333 |
2 |
| WT |
exp3 |
exp3 WT |
X2h |
34.066667 |
2 |
| PARP1KO |
exp3 |
exp3 PARP1KO |
X2h |
39.766667 |
2 |
| ALC1KO |
exp3 |
exp3 ALC1KO |
X2h |
54.466667 |
2 |
| ALC1KOPARP1KO |
exp3 |
exp3 ALC1KOPARP1KO |
X2h |
42.433333 |
2 |
| WT |
exp1 |
exp1 WT |
X4h |
29.566667 |
4 |
| PARP1KO |
exp1 |
exp1 PARP1KO |
X4h |
46.266667 |
4 |
| ALC1KO |
exp1 |
exp1 ALC1KO |
X4h |
59.766667 |
4 |
| ALC1KOPARP1KO |
exp1 |
exp1 ALC1KOPARP1KO |
X4h |
43.000000 |
4 |
| WT |
exp2 |
exp2 WT |
X4h |
29.733333 |
4 |
| PARP1KO |
exp2 |
exp2 PARP1KO |
X4h |
49.600000 |
4 |
| ALC1KO |
exp2 |
exp2 ALC1KO |
X4h |
59.833333 |
4 |
| ALC1KOPARP1KO |
exp2 |
exp2 ALC1KOPARP1KO |
X4h |
46.233333 |
4 |
| WT |
exp3 |
exp3 WT |
X4h |
35.433333 |
4 |
| PARP1KO |
exp3 |
exp3 PARP1KO |
X4h |
41.266667 |
4 |
| ALC1KO |
exp3 |
exp3 ALC1KO |
X4h |
52.100000 |
4 |
| ALC1KOPARP1KO |
exp3 |
exp3 ALC1KOPARP1KO |
X4h |
42.900000 |
4 |
| WT |
exp1 |
exp1 WT |
X6h |
21.233333 |
6 |
| PARP1KO |
exp1 |
exp1 PARP1KO |
X6h |
45.966667 |
6 |
| ALC1KO |
exp1 |
exp1 ALC1KO |
X6h |
55.666667 |
6 |
| ALC1KOPARP1KO |
exp1 |
exp1 ALC1KOPARP1KO |
X6h |
46.400000 |
6 |
| WT |
exp2 |
exp2 WT |
X6h |
23.000000 |
6 |
| PARP1KO |
exp2 |
exp2 PARP1KO |
X6h |
43.933333 |
6 |
| ALC1KO |
exp2 |
exp2 ALC1KO |
X6h |
57.933333 |
6 |
| ALC1KOPARP1KO |
exp2 |
exp2 ALC1KOPARP1KO |
X6h |
42.433333 |
6 |
| WT |
exp3 |
exp3 WT |
X6h |
20.666667 |
6 |
| PARP1KO |
exp3 |
exp3 PARP1KO |
X6h |
33.566667 |
6 |
| ALC1KO |
exp3 |
exp3 ALC1KO |
X6h |
55.133333 |
6 |
| ALC1KOPARP1KO |
exp3 |
exp3 ALC1KOPARP1KO |
X6h |
35.966667 |
6 |
Plot Data
library(ggplot2)
# quadratic
ggplot(dataset, aes(x=time, y=Outcome)) +
theme_bw() +
theme(panel.grid=element_blank(), text = element_text(size=14)) +
geom_smooth(method=lm, formula = y ~ poly(x,2), se=FALSE, aes(colour=genotype)) +
geom_point(aes(colour=genotype, shape=experiment), size=2) +
facet_grid(. ~ genotype) +
xlab(label = "Time post 2 Gy irradiation") +
scale_shape_manual(values=15:20) +
scale_color_manual(values=c("#000000","#000080","#808080","#800000"))

# cubic
ggplot(dataset, aes(x=time, y=Outcome)) +
theme_bw() +
theme(panel.grid=element_blank(), text = element_text(size=14)) +
geom_smooth(method=lm, formula = y ~ poly(x,3), se=FALSE, aes(colour=genotype)) +
geom_point(aes(colour=genotype, shape=experiment), size=2) +
facet_grid(. ~ genotype) +
xlab(label = "Time post 2 Gy irradiation") +
scale_shape_manual(values=15:20) +
scale_color_manual(values=c("#000000","#000080","#808080","#800000"))

dataplot <- aggregate(Outcome ~ genotype+treatment, data = dataset, mean)
dataplot$se <- aggregate(Outcome ~ genotype+treatment, data = dataset, FUN = function(x){ sd(x)/sqrt(length(x))})$Outcome
dataplot$treatment <- factor(gsub("X|.nM","", dataplot$treatment))
dataplot$genotype <- factor(dataplot$genotype, levels = c("WT","PARP1KO","ALC1KO","ALC1KOPARP1KO"))
dataplot$treatment <- relevel(dataplot$treatment, ref = "NT")
ggplot(dataplot, aes(x=treatment, y=Outcome, fill = genotype)) +
theme_bw() +
theme(panel.grid=element_blank(), text = element_text(size=14)) +
geom_bar(stat="identity", position=position_dodge()) +
geom_errorbar(aes(ymin=Outcome-se, ymax=Outcome+se),width=.2, position=position_dodge(.9)) +
ylab(label = paste0(expression("\u03B3"),"H2AX foci per cell")) +
xlab(label = "Time post 2 Gy irradiation") +
ylim(0, 70) +
scale_fill_manual(values=c("#000000","#800000","#000080","#808080"))

library(Cairo)
cairo_pdf("Figure3D_v1.pdf", width = 7, height = 4, family = "Arial")
ggplot(dataplot, aes(x=treatment, y=Outcome, fill = genotype)) +
theme_bw() +
theme(panel.grid.major=element_blank(), panel.grid.minor=element_blank(),
axis.line = element_line(colour = "black"), text = element_text(size=14),
panel.border = element_blank(), panel.background = element_blank()) +
geom_bar(stat="identity", position=position_dodge()) +
geom_errorbar(aes(ymin=Outcome-se, ymax=Outcome+se),width=.2, position=position_dodge(.9)) +
ylab(label = paste0(expression("\u03B3"),"H2AX foci per cell")) +
xlab(label = "Time post 2 Gy irradiation (h)") +
ylim(0, 70) +
scale_fill_manual(values=c("#000000","#800000","#000080","#808080"))
dev.off()
## quartz_off_screen
## 2
dataset$genotype <- factor(dataset$genotype, levels = c("WT","PARP1KO","ALC1KO","ALC1KOPARP1KO"))
dataset$treatment <- factor(gsub("X","", dataset$treatment))
dataset$treatment <- relevel(dataset$treatment, ref = "NT")
cairo_pdf("Figure3D_v2.pdf", width = 7, height = 4, family = "Arial")
ggplot(dataplot, aes(x=treatment, y=Outcome, fill = genotype)) +
theme_bw() +
theme(panel.grid.major=element_blank(), panel.grid.minor=element_blank(),
axis.line = element_line(colour = "black"), text = element_text(size=14),
panel.border = element_blank(), panel.background = element_blank()) +
geom_bar(stat="identity", position=position_dodge(), colour="black") +
ylab(label = paste0(expression("\u03B3"),"H2AX foci per cell")) +
xlab(label = "Time post 2 Gy irradiation (h)") +
ylim(0, 70) +
scale_fill_manual(values=c("#00000080","#80000080","#00008080","#80808080")) +
geom_jitter(data = dataset, cex=1, position=position_dodge(0.9), aes(colour=genotype)) +
scale_color_manual(values=c("#000000","#800000","#000080","#808080")) +
geom_errorbar(aes(ymin=Outcome-se, ymax=Outcome+se),width=.2, position=position_dodge(.9))
cairo_pdf("Figure3D_v3.pdf", width = 7, height = 4, family = "Arial")
ggplot(dataset, aes(x=treatment, y=Outcome, color=genotype)) +
theme_bw() +
theme(panel.grid.major=element_blank(), panel.grid.minor=element_blank(),
axis.line = element_line(colour = "black"), text = element_text(size=14),
panel.border = element_blank(), panel.background = element_blank()) +
geom_jitter(cex=1, position=position_dodge(1)) +
stat_summary(fun.data=mean_se, fun.args = list(mult=1), geom="errorbar", width=0.25, aes(colour=genotype), position=position_dodge(1)) +
stat_summary(fun.y=mean, geom="crossbar", width=0.5, aes(colour=genotype), position=position_dodge(1)) +
ylab(label = paste0(expression("\u03B3"),"H2AX foci per cell")) +
xlab(label = "Time post 2 Gy irradiation (h)") +
ylim(0, 70) +
scale_fill_manual(values=c("#00000080","#80000080","#00008080","#80808080")) +
scale_color_manual(values=c("#000000","#800000","#000080","#808080"))
dev.off()
## quartz_off_screen
## 2
Linear Model
library(MASS)
library(DHARMa)
library(lme4)
library(lmerTest)
library(bbmle)
fit1 <- lm(Outcome ~ experiment+time*genotype, data = dataset)
print(summary(fit1))
##
## Call:
## lm(formula = Outcome ~ experiment + time * genotype, data = dataset)
##
## Residuals:
## Min 1Q Median 3Q Max
## -17.090 -11.403 -2.183 11.781 25.076
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 19.66444 7.29796 2.695 0.0104 *
## experimentexp2 -0.08542 4.92028 -0.017 0.9862
## experimentexp3 -4.06458 4.92028 -0.826 0.4139
## time 1.91556 1.79663 1.066 0.2931
## genotypePARP1KO 2.63111 9.50689 0.277 0.7835
## genotypeALC1KO 3.95889 9.50689 0.416 0.6794
## genotypeALC1KOPARP1KO 2.88556 9.50689 0.304 0.7631
## time:genotypePARP1KO 3.01556 2.54082 1.187 0.2427
## time:genotypeALC1KO 5.40444 2.54082 2.127 0.0400 *
## time:genotypeALC1KOPARP1KO 2.95111 2.54082 1.161 0.2527
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 13.92 on 38 degrees of freedom
## Multiple R-squared: 0.5491, Adjusted R-squared: 0.4423
## F-statistic: 5.141 on 9 and 38 DF, p-value: 0.0001444
cat("AIC: ", AIC(fit1))
## AIC: 399.7809
simres <- simulateResiduals(fittedModel = fit1)
plot(simres)

fit2 <- lm(Outcome ~ experiment+poly(time,2)*genotype, data = dataset)
print(summary(fit2))
##
## Call:
## lm(formula = Outcome ~ experiment + poly(time, 2) * genotype,
## data = dataset)
##
## Residuals:
## Min 1Q Median 3Q Max
## -10.4611 -4.2647 -0.0672 3.4001 13.5978
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 25.41111 2.28095 11.141 6.84e-13 ***
## experimentexp2 -0.08542 2.28095 -0.037 0.9703
## experimentexp3 -4.06458 2.28095 -1.782 0.0837 .
## poly(time, 2)1 29.67566 12.90300 2.300 0.0277 *
## poly(time, 2)2 -65.08662 12.90300 -5.044 1.51e-05 ***
## genotypePARP1KO 11.67778 2.63381 4.434 9.19e-05 ***
## genotypeALC1KO 20.17222 2.63381 7.659 6.66e-09 ***
## genotypeALC1KOPARP1KO 11.73889 2.63381 4.457 8.59e-05 ***
## poly(time, 2)1:genotypePARP1KO 46.71679 18.24760 2.560 0.0151 *
## poly(time, 2)2:genotypePARP1KO -14.43376 18.24760 -0.791 0.4344
## poly(time, 2)1:genotypeALC1KO 83.72529 18.24760 4.588 5.83e-05 ***
## poly(time, 2)2:genotypeALC1KO -18.16729 18.24760 -0.996 0.3265
## poly(time, 2)1:genotypeALC1KOPARP1KO 45.71842 18.24760 2.505 0.0172 *
## poly(time, 2)2:genotypeALC1KOPARP1KO -14.20282 18.24760 -0.778 0.4418
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 6.452 on 34 degrees of freedom
## Multiple R-squared: 0.9133, Adjusted R-squared: 0.8801
## F-statistic: 27.55 on 13 and 34 DF, p-value: 2.899e-14
cat("AIC: ", AIC(fit2))
## AIC: 328.6398
simres <- simulateResiduals(fittedModel = fit2)
plot(simres)

fit3 <- lm(Outcome ~ experiment+poly(time,3)*genotype, data = dataset)
print(summary(fit3))
##
## Call:
## lm(formula = Outcome ~ experiment + poly(time, 3) * genotype,
## data = dataset)
##
## Residuals:
## Min 1Q Median 3Q Max
## -6.2076 -1.9667 0.0132 1.7253 7.1944
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 25.41111 1.21180 20.970 < 2e-16 ***
## experimentexp2 -0.08542 1.21180 -0.070 0.944274
## experimentexp3 -4.06458 1.21180 -3.354 0.002170 **
## poly(time, 3)1 29.67566 6.85499 4.329 0.000154 ***
## poly(time, 3)2 -65.08662 6.85499 -9.495 1.51e-10 ***
## poly(time, 3)3 19.41656 6.85499 2.832 0.008176 **
## genotypePARP1KO 11.67778 1.39927 8.346 2.58e-09 ***
## genotypeALC1KO 20.17222 1.39927 14.416 5.01e-15 ***
## genotypeALC1KOPARP1KO 11.73889 1.39927 8.389 2.31e-09 ***
## poly(time, 3)1:genotypePARP1KO 46.71679 9.69443 4.819 3.88e-05 ***
## poly(time, 3)2:genotypePARP1KO -14.43376 9.69443 -1.489 0.146959
## poly(time, 3)3:genotypePARP1KO 13.65011 9.69443 1.408 0.169399
## poly(time, 3)1:genotypeALC1KO 83.72529 9.69443 8.636 1.24e-09 ***
## poly(time, 3)2:genotypeALC1KO -18.16729 9.69443 -1.874 0.070701 .
## poly(time, 3)3:genotypeALC1KO 13.13372 9.69443 1.355 0.185606
## poly(time, 3)1:genotypeALC1KOPARP1KO 45.71842 9.69443 4.716 5.19e-05 ***
## poly(time, 3)2:genotypeALC1KOPARP1KO -14.20282 9.69443 -1.465 0.153310
## poly(time, 3)3:genotypeALC1KOPARP1KO 22.06740 9.69443 2.276 0.030131 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.427 on 30 degrees of freedom
## Multiple R-squared: 0.9784, Adjusted R-squared: 0.9662
## F-statistic: 79.96 on 17 and 30 DF, p-value: < 2.2e-16
cat("AIC: ", AIC(fit3))
## AIC: 269.9136
simres <- simulateResiduals(fittedModel = fit3)
plot(simres)

Compare Results
ICtab(fit1,fit2,fit3,
base=T)
## AIC dAIC df
## fit3 269.9 0.0 19
## fit2 328.6 58.7 15
## fit1 399.8 129.9 11
Final Result
fit <- fit3
output <- coef(summary(fit))
output <- output[grep("time", rownames(output)),]
rownames(output) <- gsub("poly\\(|, [1-3]\\)","", rownames(output) )
rownames(output) <- gsub("genotype", paste0(" ",levels(dataset$genotype)[1], " vs. "), rownames(output))
rownames(output)[!(grepl("vs", rownames(output)))] <- paste(rownames(output)[!(grepl("vs", rownames(output)))], levels(dataset$genotype)[1], sep = " in " )
# suggested result table
kable(output, row.names = T)
| time1 in WT |
29.67566 |
6.854994 |
4.329057 |
0.0001535 |
| time2 in WT |
-65.08662 |
6.854994 |
-9.494774 |
0.0000000 |
| time3 in WT |
19.41656 |
6.854994 |
2.832469 |
0.0081758 |
| time1: WT vs. PARP1KO |
46.71679 |
9.694426 |
4.818933 |
0.0000388 |
| time2: WT vs. PARP1KO |
-14.43376 |
9.694426 |
-1.488872 |
0.1469593 |
| time3: WT vs. PARP1KO |
13.65011 |
9.694426 |
1.408038 |
0.1693989 |
| time1: WT vs. ALC1KO |
83.72529 |
9.694426 |
8.636437 |
0.0000000 |
| time2: WT vs. ALC1KO |
-18.16729 |
9.694426 |
-1.873993 |
0.0707005 |
| time3: WT vs. ALC1KO |
13.13372 |
9.694426 |
1.354770 |
0.1856062 |
| time1: WT vs. ALC1KOPARP1KO |
45.71842 |
9.694426 |
4.715949 |
0.0000519 |
| time2: WT vs. ALC1KOPARP1KO |
-14.20282 |
9.694426 |
-1.465050 |
0.1533098 |
| time3: WT vs. ALC1KOPARP1KO |
22.06740 |
9.694426 |
2.276298 |
0.0301313 |
write.table(output, file = "Figure3D_Stats_Ref_WT.txt", quote = F, sep = "\t", row.names = T, col.names = NA)
# re-fit with ALC1KO reference
dataset$genotype <- relevel(dataset$genotype, ref = "ALC1KO")
fit <- lm(Outcome ~ experiment+poly(time,3)*genotype, data = dataset)
output <- coef(summary(fit))
output <- output[grep("time", rownames(output)),]
rownames(output) <- gsub("poly\\(|, [1-3]\\)","", rownames(output) )
rownames(output) <- gsub("genotype", paste0(" ",levels(dataset$genotype)[1], " vs. "), rownames(output))
rownames(output)[!(grepl("vs", rownames(output)))] <- paste(rownames(output)[!(grepl("vs", rownames(output)))], levels(dataset$genotype)[1], sep = " in " )
# suggested result table
kable(output, row.names = T)
| time1 in ALC1KO |
113.4009524 |
6.854994 |
16.5428227 |
0.0000000 |
| time2 in ALC1KO |
-83.2539088 |
6.854994 |
-12.1450008 |
0.0000000 |
| time3 in ALC1KO |
32.5502734 |
6.854994 |
4.7484028 |
0.0000474 |
| time1: ALC1KO vs. WT |
-83.7252933 |
9.694426 |
-8.6364367 |
0.0000000 |
| time2: ALC1KO vs. WT |
18.1672885 |
9.694426 |
1.8739933 |
0.0707005 |
| time3: ALC1KO vs. WT |
-13.1337169 |
9.694426 |
-1.3547700 |
0.1856062 |
| time1: ALC1KO vs. PARP1KO |
-37.0085075 |
9.694426 |
-3.8175039 |
0.0006289 |
| time2: ALC1KO vs. PARP1KO |
3.7335317 |
9.694426 |
0.3851215 |
0.7028645 |
| time3: ALC1KO vs. PARP1KO |
0.5163978 |
9.694426 |
0.0532675 |
0.9578719 |
| time1: ALC1KO vs. ALC1KOPARP1KO |
-38.0068766 |
9.694426 |
-3.9204877 |
0.0004749 |
| time2: ALC1KO vs. ALC1KOPARP1KO |
3.9644719 |
9.694426 |
0.4089434 |
0.6854860 |
| time3: ALC1KO vs. ALC1KOPARP1KO |
8.9336816 |
9.694426 |
0.9215277 |
0.3641269 |
write.table(output, file = "Figure3D_Stats_Ref_ALC1KO.txt", quote = F, sep = "\t", row.names = T, col.names = NA)
ANOVA
fit3a <- lm(Outcome ~ experiment+poly(time,3)*genotype, data = dataset)
fit3b <- lm(Outcome ~ experiment+poly(time,3)+genotype, data = dataset)
anova(fit3a, fit3b)
## Analysis of Variance Table
##
## Model 1: Outcome ~ experiment + poly(time, 3) * genotype
## Model 2: Outcome ~ experiment + poly(time, 3) + genotype
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 30 352.43
## 2 39 1343.95 -9 -991.51 9.3778 1.297e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
New analysis
- fit model for each time point independently
dataset$genotype <- relevel(dataset$genotype, ref = "WT")
for(i in 2:length(levels(dataset$treatment))){
dataset_sub <- dataset[dataset$treatment == levels(dataset$treatment)[i],]
for(g in seq_along(levels(dataset_sub$genotype))){
cat("Time: ", levels(dataset$treatment)[i], "\n",
"Genotype: ", levels(dataset$genotype)[g], "\n")
dataset_sub$genotype <- relevel(dataset_sub$genotype, ref = levels(dataset_sub$genotype)[g])
fit_sub <- lm(Outcome ~ genotype, data = dataset_sub)
print(summary(fit_sub))
cat("AIC: ", AIC(fit_sub))
simres <- simulateResiduals(fittedModel = fit_sub)
plot(simres)
if(i == 2 & g == 1){
output <- coef(summary(fit_sub))
output <- output[grep("genotype", rownames(output)),]
rownames(output) <- paste0(gsub("genotype", paste0(" ",levels(dataset$genotype)[g], " vs. "), rownames(output)),
" in Time ", levels(dataset$treatment)[i])
} else {
outtmp <- coef(summary(fit_sub))
outtmp <- outtmp[grep("genotype", rownames(outtmp)),]
rownames(outtmp) <- paste0(gsub("genotype", paste0(" ",levels(dataset$genotype)[g], " vs. "), rownames(outtmp)),
" in Time ", levels(dataset$treatment)[i])
output <- rbind(output,outtmp)
if(g == length(levels(dataset_sub$genotype)) & i < length(levels(dataset$treatment))){
output <- rbind(output, " ", colnames(output))
}
}
}
}
## Time: 2h
## Genotype: WT
##
## Call:
## lm(formula = Outcome ~ genotype, data = dataset_sub)
##
## Residuals:
## Min 1Q Median 3Q Max
## -8.8889 -1.0750 0.3222 2.0889 8.5778
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 35.267 3.266 10.797 4.78e-06 ***
## genotypeALC1KO 19.933 4.619 4.315 0.00256 **
## genotypePARP1KO 13.389 4.619 2.898 0.01994 *
## genotypeALC1KOPARP1KO 15.111 4.619 3.271 0.01134 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 5.658 on 8 degrees of freedom
## Multiple R-squared: 0.7191, Adjusted R-squared: 0.6137
## F-statistic: 6.826 on 3 and 8 DF, p-value: 0.01349
##
## AIC: 80.78097

## Time: 2h
## Genotype: ALC1KO
##
## Call:
## lm(formula = Outcome ~ genotype, data = dataset_sub)
##
## Residuals:
## Min 1Q Median 3Q Max
## -8.8889 -1.0750 0.3222 2.0889 8.5778
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 55.200 3.266 16.899 1.52e-07 ***
## genotypeWT -19.933 4.619 -4.315 0.00256 **
## genotypePARP1KO -6.544 4.619 -1.417 0.19431
## genotypeALC1KOPARP1KO -4.822 4.619 -1.044 0.32705
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 5.658 on 8 degrees of freedom
## Multiple R-squared: 0.7191, Adjusted R-squared: 0.6137
## F-statistic: 6.826 on 3 and 8 DF, p-value: 0.01349
##
## AIC: 80.78097

## Time: 2h
## Genotype: PARP1KO
##
## Call:
## lm(formula = Outcome ~ genotype, data = dataset_sub)
##
## Residuals:
## Min 1Q Median 3Q Max
## -8.8889 -1.0750 0.3222 2.0889 8.5778
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 48.656 3.266 14.896 4.07e-07 ***
## genotypeALC1KO 6.544 4.619 1.417 0.1943
## genotypeWT -13.389 4.619 -2.898 0.0199 *
## genotypeALC1KOPARP1KO 1.722 4.619 0.373 0.7190
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 5.658 on 8 degrees of freedom
## Multiple R-squared: 0.7191, Adjusted R-squared: 0.6137
## F-statistic: 6.826 on 3 and 8 DF, p-value: 0.01349
##
## AIC: 80.78097

## Time: 2h
## Genotype: ALC1KOPARP1KO
##
## Call:
## lm(formula = Outcome ~ genotype, data = dataset_sub)
##
## Residuals:
## Min 1Q Median 3Q Max
## -8.8889 -1.0750 0.3222 2.0889 8.5778
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 50.378 3.266 15.423 3.11e-07 ***
## genotypePARP1KO -1.722 4.619 -0.373 0.7190
## genotypeALC1KO 4.822 4.619 1.044 0.3270
## genotypeWT -15.111 4.619 -3.271 0.0113 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 5.658 on 8 degrees of freedom
## Multiple R-squared: 0.7191, Adjusted R-squared: 0.6137
## F-statistic: 6.826 on 3 and 8 DF, p-value: 0.01349
##
## AIC: 80.78097

## Time: 4h
## Genotype: WT
##
## Call:
## lm(formula = Outcome ~ genotype, data = dataset_sub)
##
## Residuals:
## Min 1Q Median 3Q Max
## -5.1333 -1.8861 -0.2444 2.5500 3.8889
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 31.578 2.084 15.154 3.56e-07 ***
## genotypeALC1KO 25.656 2.947 8.706 2.36e-05 ***
## genotypePARP1KO 14.133 2.947 4.796 0.00136 **
## genotypeALC1KOPARP1KO 12.467 2.947 4.230 0.00288 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.609 on 8 degrees of freedom
## Multiple R-squared: 0.9049, Adjusted R-squared: 0.8693
## F-statistic: 25.39 on 3 and 8 DF, p-value: 0.0001931
##
## AIC: 69.99344

## Time: 4h
## Genotype: ALC1KO
##
## Call:
## lm(formula = Outcome ~ genotype, data = dataset_sub)
##
## Residuals:
## Min 1Q Median 3Q Max
## -5.1333 -1.8861 -0.2444 2.5500 3.8889
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 57.233 2.084 27.465 3.33e-09 ***
## genotypeWT -25.656 2.947 -8.706 2.36e-05 ***
## genotypePARP1KO -11.522 2.947 -3.910 0.00448 **
## genotypeALC1KOPARP1KO -13.189 2.947 -4.475 0.00207 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.609 on 8 degrees of freedom
## Multiple R-squared: 0.9049, Adjusted R-squared: 0.8693
## F-statistic: 25.39 on 3 and 8 DF, p-value: 0.0001931
##
## AIC: 69.99344

## Time: 4h
## Genotype: PARP1KO
##
## Call:
## lm(formula = Outcome ~ genotype, data = dataset_sub)
##
## Residuals:
## Min 1Q Median 3Q Max
## -5.1333 -1.8861 -0.2444 2.5500 3.8889
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 45.711 2.084 21.936 1.97e-08 ***
## genotypeALC1KO 11.522 2.947 3.910 0.00448 **
## genotypeWT -14.133 2.947 -4.796 0.00136 **
## genotypeALC1KOPARP1KO -1.667 2.947 -0.566 0.58721
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.609 on 8 degrees of freedom
## Multiple R-squared: 0.9049, Adjusted R-squared: 0.8693
## F-statistic: 25.39 on 3 and 8 DF, p-value: 0.0001931
##
## AIC: 69.99344

## Time: 4h
## Genotype: ALC1KOPARP1KO
##
## Call:
## lm(formula = Outcome ~ genotype, data = dataset_sub)
##
## Residuals:
## Min 1Q Median 3Q Max
## -5.1333 -1.8861 -0.2444 2.5500 3.8889
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 44.044 2.084 21.136 2.64e-08 ***
## genotypePARP1KO 1.667 2.947 0.566 0.58721
## genotypeALC1KO 13.189 2.947 4.475 0.00207 **
## genotypeWT -12.467 2.947 -4.230 0.00288 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.609 on 8 degrees of freedom
## Multiple R-squared: 0.9049, Adjusted R-squared: 0.8693
## F-statistic: 25.39 on 3 and 8 DF, p-value: 0.0001931
##
## AIC: 69.99344

## Time: 6h
## Genotype: WT
##
## Call:
## lm(formula = Outcome ~ genotype, data = dataset_sub)
##
## Residuals:
## Min 1Q Median 3Q Max
## -7.5889 -1.0028 0.2167 1.9611 4.8111
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 21.633 2.511 8.616 2.55e-05 ***
## genotypeALC1KO 34.611 3.551 9.747 1.03e-05 ***
## genotypePARP1KO 19.522 3.551 5.498 0.000575 ***
## genotypeALC1KOPARP1KO 19.967 3.551 5.623 0.000497 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.349 on 8 degrees of freedom
## Multiple R-squared: 0.9231, Adjusted R-squared: 0.8942
## F-statistic: 31.99 on 3 and 8 DF, p-value: 8.358e-05
##
## AIC: 74.46729

## Time: 6h
## Genotype: ALC1KO
##
## Call:
## lm(formula = Outcome ~ genotype, data = dataset_sub)
##
## Residuals:
## Min 1Q Median 3Q Max
## -7.5889 -1.0028 0.2167 1.9611 4.8111
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 56.244 2.511 22.400 1.67e-08 ***
## genotypeWT -34.611 3.551 -9.747 1.03e-05 ***
## genotypePARP1KO -15.089 3.551 -4.249 0.00280 **
## genotypeALC1KOPARP1KO -14.644 3.551 -4.124 0.00333 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.349 on 8 degrees of freedom
## Multiple R-squared: 0.9231, Adjusted R-squared: 0.8942
## F-statistic: 31.99 on 3 and 8 DF, p-value: 8.358e-05
##
## AIC: 74.46729

## Time: 6h
## Genotype: PARP1KO
##
## Call:
## lm(formula = Outcome ~ genotype, data = dataset_sub)
##
## Residuals:
## Min 1Q Median 3Q Max
## -7.5889 -1.0028 0.2167 1.9611 4.8111
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 41.1556 2.5109 16.391 1.93e-07 ***
## genotypeALC1KO 15.0889 3.5509 4.249 0.002802 **
## genotypeWT -19.5222 3.5509 -5.498 0.000575 ***
## genotypeALC1KOPARP1KO 0.4444 3.5509 0.125 0.903482
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.349 on 8 degrees of freedom
## Multiple R-squared: 0.9231, Adjusted R-squared: 0.8942
## F-statistic: 31.99 on 3 and 8 DF, p-value: 8.358e-05
##
## AIC: 74.46729

## Time: 6h
## Genotype: ALC1KOPARP1KO
##
## Call:
## lm(formula = Outcome ~ genotype, data = dataset_sub)
##
## Residuals:
## Min 1Q Median 3Q Max
## -7.5889 -1.0028 0.2167 1.9611 4.8111
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 41.6000 2.5109 16.568 1.78e-07 ***
## genotypePARP1KO -0.4444 3.5509 -0.125 0.903482
## genotypeALC1KO 14.6444 3.5509 4.124 0.003325 **
## genotypeWT -19.9667 3.5509 -5.623 0.000497 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.349 on 8 degrees of freedom
## Multiple R-squared: 0.9231, Adjusted R-squared: 0.8942
## F-statistic: 31.99 on 3 and 8 DF, p-value: 8.358e-05
##
## AIC: 74.46729

# suggested result table
kable(output, row.names = T)
| WT vs. ALC1KO in Time 2h |
19.9333333366666 |
4.61941688316092 |
4.31511895133979 |
0.00256288937859708 |
| WT vs. PARP1KO in Time 2h |
13.38888889 |
4.61941688316092 |
2.89839372125219 |
0.0199407600005576 |
| WT vs. ALC1KOPARP1KO in Time 2h |
15.11111111 |
4.61941688316092 |
3.27121614961496 |
0.0113355916961444 |
| ALC1KO vs. WT in Time 2h |
-19.9333333366666 |
4.61941688316091 |
-4.31511895133979 |
0.00256288937859707 |
| ALC1KO vs. PARP1KO in Time 2h |
-6.54444444666667 |
4.61941688316091 |
-1.41672523008759 |
0.194305008686908 |
| ALC1KO vs. ALC1KOPARP1KO in Time 2h |
-4.82222222666667 |
4.61941688316091 |
-1.04390280172483 |
0.32704543157344 |
| PARP1KO vs. ALC1KO in Time 2h |
6.54444444666666 |
4.61941688316092 |
1.41672523008759 |
0.194305008686909 |
| PARP1KO vs. WT in Time 2h |
-13.38888889 |
4.61941688316092 |
-2.89839372125219 |
0.0199407600005575 |
| PARP1KO vs. ALC1KOPARP1KO in Time 2h |
1.72222222 |
4.61941688316092 |
0.372822428362763 |
0.718964761240989 |
| ALC1KOPARP1KO vs. PARP1KO in Time 2h |
-1.72222222 |
4.61941688316091 |
-0.372822428362763 |
0.71896476124099 |
| ALC1KOPARP1KO vs. ALC1KO in Time 2h |
4.82222222666667 |
4.61941688316092 |
1.04390280172483 |
0.32704543157344 |
| ALC1KOPARP1KO vs. WT in Time 2h |
-15.11111111 |
4.61941688316091 |
-3.27121614961496 |
0.0113355916961444 |
|
|
|
|
|
|
Estimate |
Std. Error |
t value |
Pr(>|t|) |
| WT vs. ALC1KO in Time 4h |
25.6555555566666 |
2.94700098560667 |
8.70564878735025 |
2.36474253843872e-05 |
| WT vs. PARP1KO in Time 4h |
14.1333333366667 |
2.94700098560667 |
4.79583597212716 |
0.00136277121174079 |
| WT vs. ALC1KOPARP1KO in Time 4h |
12.4666666666667 |
2.94700098560667 |
4.23028927630314 |
0.00287522676347967 |
| ALC1KO vs. WT in Time 4h |
-25.6555555566667 |
2.94700098560666 |
-8.70564878735025 |
2.3647425384387e-05 |
| ALC1KO vs. PARP1KO in Time 4h |
-11.52222222 |
2.94700098560666 |
-3.90981281522309 |
0.0044822300200186 |
| ALC1KO vs. ALC1KOPARP1KO in Time 4h |
-13.18888889 |
2.94700098560666 |
-4.47535951104711 |
0.00206847786110044 |
| PARP1KO vs. ALC1KO in Time 4h |
11.52222222 |
2.94700098560667 |
3.90981281522308 |
0.00448223002001862 |
| PARP1KO vs. WT in Time 4h |
-14.1333333366667 |
2.94700098560667 |
-4.79583597212717 |
0.00136277121174079 |
| PARP1KO vs. ALC1KOPARP1KO in Time 4h |
-1.66666667 |
2.94700098560667 |
-0.565546695824028 |
0.587208761244944 |
| ALC1KOPARP1KO vs. PARP1KO in Time 4h |
1.66666667 |
2.94700098560666 |
0.565546695824028 |
0.587208761244943 |
| ALC1KOPARP1KO vs. ALC1KO in Time 4h |
13.18888889 |
2.94700098560667 |
4.47535951104711 |
0.00206847786110044 |
| ALC1KOPARP1KO vs. WT in Time 4h |
-12.4666666666667 |
2.94700098560666 |
-4.23028927630314 |
0.00287522676347967 |
|
|
|
|
|
|
Estimate |
Std. Error |
t value |
Pr(>|t|) |
| WT vs. ALC1KO in Time 6h |
34.61111111 |
3.55089103646649 |
9.747162262811 |
1.02732215316808e-05 |
| WT vs. PARP1KO in Time 6h |
19.5222222233333 |
3.55089103646649 |
5.49783759142325 |
0.000575210034763358 |
| WT vs. ALC1KOPARP1KO in Time 6h |
19.9666666666667 |
3.55089103646649 |
5.62300179352606 |
0.000496757913056301 |
| ALC1KO vs. WT in Time 6h |
-34.61111111 |
3.55089103646648 |
-9.74716226281101 |
1.02732215316807e-05 |
| ALC1KO vs. PARP1KO in Time 6h |
-15.0888888866667 |
3.55089103646648 |
-4.24932467138775 |
0.00280172251246753 |
| ALC1KO vs. ALC1KOPARP1KO in Time 6h |
-14.6444444433333 |
3.55089103646648 |
-4.12416046928495 |
0.00332508800054009 |
| PARP1KO vs. ALC1KO in Time 6h |
15.0888888866667 |
3.55089103646649 |
4.24932467138775 |
0.00280172251246755 |
| PARP1KO vs. WT in Time 6h |
-19.5222222233333 |
3.55089103646649 |
-5.49783759142326 |
0.000575210034763353 |
| PARP1KO vs. ALC1KOPARP1KO in Time 6h |
0.44444444333333 |
3.55089103646649 |
0.125164202102805 |
0.903481679641302 |
| ALC1KOPARP1KO vs. PARP1KO in Time 6h |
-0.444444443333329 |
3.55089103646648 |
-0.125164202102805 |
0.903481679641302 |
| ALC1KOPARP1KO vs. ALC1KO in Time 6h |
14.6444444433333 |
3.55089103646649 |
4.12416046928495 |
0.0033250880005401 |
| ALC1KOPARP1KO vs. WT in Time 6h |
-19.9666666666667 |
3.55089103646648 |
-5.62300179352606 |
0.000496757913056298 |
write.table(output, file = "Figure3D_Stats_New_All.txt", quote = F, sep = "\t", row.names = T, col.names = NA)